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,^ ■ Abstract 

^^ I For nonlinear reduced-order models, especially for those with high-order polynomial non- 

^ ■ linearities or non-polynomial nonlinearities, the computational complexity still depends on the 

4— i I dimension of the original dynamical system. As a result, the reduced-order model loses its 

2 ■ computational efficiency, which, however, is its the most significant advantage. Nonlinear di- 

G I mensional reduction methods, such as the discrete empirical interpolation method, have been 

widely used to evaluate the nonlinear terms at a low cost. But when the finite element method 
^S| I is utilized for the spatial discretization, nonlinear snapshot generation requires inner products 

^ . to be fulfilled, which costs lots of off-line time. Numerical integrations are also needed over ele- 

wN ■ ments sharing the selected interpolation points during the simulation, which keeps on-line time 

^>l . high. To overcome these issues and develop an efficient finite element discretization algorithm, 

^^ I in this paper, we extend the finite element method with interpolated coefficients, also known as 

the group finite element method or the product approximation, to nonlinear reduced-order mod- 
els. The proposed approach approximates the nonlinear function in the reduced-order model 
^^ . by its finite element interpolation, which makes coefficient matrices of the nonlinear terms pre- 

^-H ' computable and, thus, leads to great savings in the computational efforts. Due to the separation 

K^ . of spatial and temporal variables in the finite element interpolation, the discrete empirical in- 

• i-H ' terpolation method can be directly applied on the nonlinear functions in the same manner as 

/^ . that in the finite difference setting. Therefore, the main computational hurdles when applying 

the discrete empirical interpolation method in the finite element context are conquered. We 
also establish a rigorous asymptotic error estimation, which shows that the proposed approach 
achieves the same accuracy as that of the standard finite element method under certain smooth- 
ness assumptions of the nonlinear functions. Several numerical tests are presented to validate 
the proposed method and verify the theoretical results. 

Keywords: Nonlinear model reduction, finite element method with interpolated coefficients, 
proper orthogonal decomposition, discrete empirical interpolation method 

1 Introduction 

Control and optimization problems in realistic engineering applications often require repeated nu- 
merical simulations of large-scale dynamical systems. If a fast or real-time control strategy is 



o 



U 



*Eniail: wangzhu@ima.umn.edu. URL: http: //www. ima.uinn.edu/~wangzliu/ 



Nonlinear Model Reduction for Semilinear Parabolic Equations Z. Wang 



desired, a brute force direct numerical simulation (DNS) is impractical. Therefore, the proper 
orthogonal decomposition method (POD) has been frequently used to generate a reduced-order 
model (ROM) that can be utilized as an alternative of the original dynamical system (hereinafter 
referred to as the full-order model) [32] . Such a ROM only contains a handful of degrees of free- 
dom (DOF), yet is computationally feasible and free of storage issues. The POD method has 
been successfully applied in many scientific and engineering problems (see, e.g., a brief summa- 
rization in [71]). However, for complex dynamical systems, the original promise of the POD as 
an efficient yet accurate approximation remains to be fulfilled. On the one hand, it may achieve 
erroneous results even when POD basis functions capture most of the system energy [4]. On the 
other hand, the efficiency is severely limited by the nonlinearity of the system [18]. For treating the 
first issue, research has been done in two main directions: (i) to improve the POD basis functions 
[1, 6, 7, 9, 12, 13, 14, 24, 38, 45, 51]; (ii) to improve the ROM by modeling the effect of discarded 
POD basis [8, 10, 22, 23, 27, 29, 30, 36, 37, 40, 48, 52, 53, 54, 55, 56, 57, 58, 59, 60, 65, 70, 71, 72, 73]. 
To our knowledge, the research for seeking a reliable reduced-order modeling approach for general 
complex systems is still active, however, beyond the scopes of this paper. Instead, we will focus on 
the second issue and develop an efficient algorithm for the nonlinear ROM. 

Indeed, the computational efficiency of POD-ROMs relies on two key components: (i) the 
dimension of the ROM is extremely small; (ii) the vectors and matrices in the reduced system can 
be precomputed. For linear systems or nonlinear systems with low order polynomial nonlinearities, 
both ingredients are satisfied. However, for highly nonlinear dynamical systems, the second one 
does not hold any more because matrices associated with the nonlinear terms have to be evaluated 
and assembled at each time step or iteration. Since POD basis functions are global, the evaluation 
would depend on the dimension of the full-order model. The reassembling process would greatly 
increase the computational cost. Therefore, several methods have been proposed to resolve this 
issue. 

The two-level algorithms proposed in [71] are motivated by the observation that only a small 
number of leading POD basis functions are kept in the ROM, and they have larger length scales 
than the discarded ones. If a computation on a fine mesh is employed to obtain all the POD basis, 
one should be able to use a much coarser mesh to represent the leading POD basis. Therefore, 
in two-level algorithms, nonlinear closure terms are computed on the coarse mesh. This way 
can decrease the computational cost by an order of magnitude, while achieving the same level of 
accuracy as the simulation on the fine mesh. However, the optimal choice of the coarse mesh still 
needs to be investigated. 

Much work devotes to approximate nonlinear terms at a few selected spatial points or within 
the neighborhoods around the points. The trajectory piecewise-linear method (TPWL), presented 
in [61, 62], reduces a nonlinear model to a weighted sum of linearized models at selected points 
along a state trajectory. The missing points estimation method (MPE) was developed in [2, 3, 74]. 
In that approach, the full-order system was first reduced by choosing equations only associated 
with selected spatial points, restricting the POD basis onto these points, and then projecting the 
extracted system onto the space spanned by the POD basis. The empirical interpolation method 
(EIM) was first proposed in [5] for approximating non-affine parameter dependence functions to 
enable an efficient offiine-online computational strategy, and then was further applied to approx- 
imate nonlinear functions in [31]. The method selects interpolation points by greedy algorithms 
guided by a posteriori error estimates. However, in certain problems, even the most optimal basis 
set is of large size [26], which reduces the efficiency of the algorithm. Thus, several improved 
greedy algorithms have been proposed in [33]. The best-points interpolation method (BPIM) was 
introduced in [49], which determines interpolation points from a least-square minimization prob- 
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lem. The discrete empirical interpolation method (DEIM) was introduced in [18] and analyzed in 
[17]. The method combines projection with interpolation and chooses interpolation points from the 
POD basis of the nonlinear functions. Since a certain coefficient matrix can be precomputed when 
approximating nonlinear functions, the complexity of the POD-ROM reduces to be proportional 
to the number of selected spatial indices. 

Another method is the group POD method (GPOD) proposed in [25], which extends the group 
finite element method to the POD setting. In that paper, the authors considered dynamical systems 
with polynomial nonlinearities, such as the Burgers equations. The POD approximation of the 
quadratic nonlinear term was rewritten in a group format. If r POD basis functions are used in the 
ROM, the GPOD requires r^ — ^r^ — ^r fiops less than the standard POD implementation in the 
computation of the nonlinear term. However, the approach is limited to polynomial nonlinearities. 

In this report, we focus on developing an efficient finite element (FE) discretization algorithm 
for nonlinear ROMs. In particular, we are interested in applying the DEIM to reduce the intensive 
computational efforts for evaluating the nonlinear terms. However, in the FE setting, there are 
two major issues that degrade the effectiveness of the DEIM: (i) generating nonlinear snapshots, 
which are to be used for seeking the nonlinear POD basis, requires calculations of the inner product 
in the nonlinear terms, which costs lots of off-line computation time; (ii) the on-line simulation 
needs evaluations of the inner product over the elements sharing selected DEIM points. Repeated 
numerical integrations will increase the on-line simulation time, especially, in cases such as complex 
ffow simulations when many DEIM points are required to achieve a good approximation. 

To overcome these hurdles, we develop the finite element method with interpolated coefficients 
(FEIC) for nonlinear POD-ROMs. This method is also known as the group finite element method 
or the product approximation, which has been successfully applied to find numerical solutions to 
nonlinear partial differential equations [19, 20, 21, 28, 39, 46, 47, 63, 69, 75]. Indeed, we replace the 
nonlinear function in the POD-ROM with its FE interpolation directly. This simple change would 
lead to great savings in the computation: First, the coefficient matrices of the nonlinear terms can 
be computed beforehand, thus, the evaluation of the nonlinear terms does not involve any numerical 
quadratures during the on-line simulation; Second, for the ROMs discretized by the new approach, 
the DEIM can be directly applied in the same manner as that in the finite difference context. In 
fact, the nonlinear snapshots become a collection of vectors of nonlinear function values, therefore, 
neither does the nonlinear data generation require any numerical integrations. The accuracy of this 
approach is, of course, restricted to the smoothness of the nonlinear functions. However, when the 
nonlinear functions possess certain smoothness ((H2) and (H3) in Section 4), the FEIC can achieve 
the same accuracy as the standard FE discretization of the ROM. Therefore, the advantages of 
the new approach over the standard FE discretization of nonlinear POD-ROMs cannot be over- 
emphasized: (i) the FEIC is easier to implement and computationally more efficient; (ii) the FEIC 
achieves the same accuracy when nonlinear functions satisfy smoothness assumptions; and (iii) 
the FEIC is more suitable for combining with the DEIM and further reduce the computational 
complexity. 

Note that the GPOD method is based on a similar idea. Distinguishing from it, the proposed 
method in this paper doesn't group any variables in terms of the POD basis. Therefore, it fits well 
for ROMs with polynomial or non-polynomial nonlinearities. The rest of this paper is organized as 
follows: a brief introduction to the POD method is presented in Section 2; the FEIC of POD-ROMs 
for semilinear parabolic equations is proposed in Section 3; a rigorous asymptotic error estimate 
is developed in Section 4; several examples are tested in Section 5 to numerically demonstrate 
the accuracy and efficiency of the new approach. Wherein, the first two examples are used for 
validation and the third example is for verification. The combination of the proposed approach 
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with the DEIM is discussed in Section 6, whose effectiveness is then illustrated by revisiting the 
first two examples used in Section 5. A few conclusions are drawn in the last section with several 
ongoing research directions we are pursuing listed. 

2 The POD Method 

In this section, we briefly introduce the (time-continuous) proper orthogonal decomposition method. 
For detailed discussions, the reader is referred to [15, 35, 43, 44, 64, 66]. 

Let T-Lhe di real Hilbert space endowed with inner product (•, •)y^ and norm || • H^^. Assume the 
data V (so-called snapshots), which is a collection of time-varying functions y{x,t) G L^(0,T;H), 
the POD method seeks a low-dimensional basis, (/pi(x), . . . , (fr{x) G 7^, that optimally approximates 
the data. Mathematically speaking, for any positive r, the POD basis is determined by minimizing 
the error between the data and its projection onto the basis, that is. 



mm 



^=1 J^ ^=1 



^dt, (1) 



subject to the conditions that {(fi,(fj)n = 5ij, I ^ i,j ^ r, where 6ij is the Kronecker delta. In 
order to solve (1), one can consider the eigenvalue problem 

Kvj^XjVj, (2) 

where if is a compact linear operator that satisfies Kvj(s) = J^ (y(-^t)^y(-^s))nVj{t) dt^ Vj is the 
j-th eigenvector, and Aj < . . . < A2 < Ai are positive eigenvalues. 
It can then be shown that the solution of (1) is given by 

^j{-) = ^ [ Vj{t)y{-,t)dt, l<j<r. (3) 

VA7 ^0 



Proposition 2.1 ([43]) Let the POD basis given by (3) be of rank r, the POD projection error 

satisfies 

rT ' 2 

y{-,t)-J2iyi-^t),^j{-))^^ji-) dt^J2^j- W 



Jo 






Remark 2.1 In practice, discrete data is always considered. The POD basis for an ensemble of 
snapshots, V = span {y{x,ti), . . . ,y{x,tM)}, is to minimize the projection error 



M „ T 2 






The solution can be obtained by solving the eigenvalue problem K Vj = Xj Vj first, where K G R^^^ 
is the snapshot correlation matrix with K^j^ = {y(-,t£),y{-,tk))^ and 1 < i,k < M, then the POD 
basis is given by (fj{-) = -^ E^i(^j)^ I/(^^^); I <j <r. 
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Remark 2.2 Finite element solutions are commonly utilized as the snapshots, that is, the data 

ridof 

y{x,t) = Yl ^i{t)hi{x), where Y(t) is the approximation solution vector at time t, h,^{x) is the 

L-th FE nodal basis, and ridof is the number of DOF in the spatial discretization. Then the POD 
basis can be written as 

ndof 

^j(^) = ^^Qi^jhiix), (6) 

L=l 

where the coefficient matrix Q is with the entry Q^j = — ^ Jq Vj{t)Yf^{t) dt for time- continuous 
data or Q^j = —j= X^^^i('^j)^^^(^^) for time-discrete data. For a detailed discussion on imple- 
menting the POD method in the FE setting, readers are referred to [42]. 

3 The POD-ROM and Finite Element Discretizations 

Let ft he a, convex domain in R^ with the smooth boundary dft, d = 1,2,3. Also let 7/i be a 
collection of quasi uniform elements that partitions the domain. The elements are line segments if 
d = 1, triangles if d = 2, and polyhedra if d = 3. The parameter h is the maximal diameter of the 
elements. Denoted by X the space L^(ri) equipped with the inner product (•,•) and norm || • ||; 
by V the Sobolev space Hq{Q) = {^'l^' G H^{^),v\qq = 0} with H^ semi-norm | • |i and H^ norm 
II • II i; and by V^ the space of piecewise continuous functions on ft that reduce to polynomials of 
degree < tti on each element of Th, which satisfies V^ C V. We assume the semilinear problems 
that we will consider in this section admit a unique solution u = u(x^ t) on the time interval [0, T], 
which ranges in V . 

We consider the equivalent variational problems of semilinear parabolic equations with homo- 
geneous boundary conditions: To find u{x^t) G V^ such that, either 

^,v] +aiu,v) + iNiu),v) = if,v), yv € V, (7) 

or 

^,v) +aiu,v) + iNiu),Vv) = if,v), ^v G V, (8) 

with the initial condition 

u{x,0) = uq{x), Vx g ri. 

Where / = f{x,t) is a source term independent with u, the bilinear form a(-, •) : "1/ x y ^ M is 
continuous and coercive, that is, there exist constants a and /3 such that 

\a(u^v)\ < a\\u\\i\\v\\i^ \/u^v^V^ (9) 

\a{u,u)\ > /3\\u\\l, \fueV. (10) 

N{u) is a nonlinear function of u. In particular, we are interested in cases in which N{u) possesses 
either a non-polynomial nonlinearity or a high order polynomial nonlinearity. Due to the similarity 
between (7) and (8), to shorten the presentation, we only discuss (7) in the sequel, but will comment 
a theoretical result of (8) in Remark 4.2 and test a problem governed by (8) in Section 5. 

Many well-established methods can be used for seeking a numerical solution to the equation 
(7). However, when repeated numerical simulations are required, direct simulations result in a 
huge computational cost and become infeasible. Therefore, the POD method has been widely used 
for generating a reduced-order model. 
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3.1 The POD-ROM 

Given snapshots consisting of either numerical solutions or experimental data, the POD basis 
functions {(fi{x), . . . , (pr{x)} are determined by (l)-(3), where T-L = L^ in our considerations. The 
associated POD approximation of u{x,t) is given by 

r 
Ur{x,t) = ^Lpj{x)aj{t), (11) 

where {aj{t)Y-Y are the sought time- varying POD basis coefficient functions. Substituting the 
POD approximation (11) into (7), applying the Galerkin method, and considering the POD basis 
functions are orthonormal, we obtain the Galerkin projection-based POD-ROM (POD-G-ROM): 
to find Ur{x^t) ^Vr = span{(/Pi, . . . , c/p^}, such that, 

"QT^^rj +a{Ur,Vr) + {N {Ur) , Vr) = (/, ^^r) , V i;^ G K, (12) 

and 

Let a(t) = [ai(t), . . . , a^(t)]T, the POD-G-ROM can be rewritten in terms of POD basis functions 
as: 

a = A + Ba + C(a), (13) 

with, for example, 

a(0) = {uo{x),^{x)). 



Where Aa; = (/,(/^/c), ^jk = -a{(fj,(fk), i^i^))k = " (^(Ej=i ^j^j W). ^/cj, for /c = 1, . . . , r, and 
$ = [(fi, . . . , (frV- The resulted dynamic system (13) is of dimension r, which is much smaller than 
the number of DOF in the full-order model. Once a is obtained, the POD approximation solution 
Ur can be recovered by (11). 

The most significant advantage of the ROM is its computational efficiency. Indeed, the matrix 
B can be precomputed. In certain cases, we can also compute matrices in C beforehand. For exam- 
ple, in Navier-Stokes equations, one can write (C(a))^ = a^C^a, where [C^]ij — — ^l=i ^^j=i{^i ' 
V(fj^(fk)' The matrices only need to be computed once and can be used repeatedly in on-line 
simulations. However, this attractive property does not hold when the nonlinear function N{u) is 
with a higher oder polynomial nonlinearity or a non-polynomial one. As a result, the computa- 
tional efficiency of the ROM decays. Especially, when the FE is used for a spatial discretization. 
Therefore, in the rest of this paper, we restrict ourselves to the FE methods of the nonlinear ROMs 
(12) and develop a new efficient FE discretization algorithm. 

3.2 Finite Element Discretizations 

Let V^ = span |(/;^, . . . , c/p^}, where (f^ is the finite element discretization of (/?j, j = 1, . . . , r. The 
standard finite element discretization of the POD-G-ROM (12) (POD-FEM) is to find u^(x,t) G 
V^, such that. 



~dt 



, v^\ + a [4, v^) + {n{u^,),v^) = (/, v^,), yv^ G K^ (14) 
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and 

Represented by Mfe the nonlinear term in (14), whose fc-th entry is as follows. 

{NFE)k = (ivfe(x,t)),^^(x) 



r ^dof ridof 

^(E E Q.Aix)ajit)), Yl Qskhsix)). (15) 

j=l L=l S=l 

Obviously, the evaluation of (15) requires numerical quadratures at each time step of the numerical 
simulation. Suppose the complexity for evaluating the nonlinear function N{u) with 9 components 
is 0{g(9)) and riq quadrature points are used in each integral evaluation, the total complexity in 
assembling AfrE at each time step (or iteration) is around O {2rnq[g{2rndof) + ^do/]) flops. 

To improve the efliciency of nonlinear ROMs, we propose a method to use the FEIC for a spatial 
discretization of the POD-G-ROM (12). The FEIC, also known as the product approximation 
technique [21] or group finite element method [28], has been used as an alternative tool of the FE 
method for solving nonlinear elliptic problems [47, 63], nonlinear parabolic problems [19, 20, 39, 46], 
and nonlinear hyperbolic problems [69, 75]. This approach replaces the nonlinear function by its 
interpolant in the finite dimensional space, which leads to great savings in the computational efforts 
while keeping the accuracy. To our knowledge, this is the first time that the FEIC is applied in 
the POD setting with a rigorous error estimate provided. 

Define the interpolation operator I^ : C{ft) -^ Sh- The interpolant of N{-) in the FE space 
satisfies 

{I^N){u{x,,t)) = N{u{x,,t)), (16) 

where Xi is a node in the finite element mesh, i = 1, . . . ^n^of- The finite element method with 
interpolated coefficients of the POD-G-ROM (12) (POD-FEIC) is the following: to find u^{x,t) G 
V^^ such that. 






, i;M + a (ix^, v^^ + (x^N{u^,),v^^ = (/, v^^ , Vi;^ G V^ , (17) 



and 



i/^(-, 0) = 4,r(^) ^ K^. Vx G n. (18) 

Different from the standard FE discretization, the nonlinear function N{u^) in (14) is replaced 
by the interpolation X^N{u^) in (17). The fc-th row of the nonlinear term in the new numerical 
discretization, Mfeic^ reads: 

{MFEic)k = ({X^N){u^Ax,t)l^l{x] 

^dof r ^dof 

^N{^Cl,jaj{t))h,{x), ^^,kn,[x, 

L=l j = l S = l 

ndof ^dof 



j=l s=l 

^dof r 

(5^/i,(x), 5^Q.fc/i.(x))iv(5^Q,,a,(t) 



L=l S = l j = l 

which can then be rewritten as follows. 



MFEic^Q^M''N{Qsi{t)), (19) 
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where M^ is the FE mass matrix with [M^] = (/i^, /i^). 

Compare with (15), the advantage of the new approach is clear: the matrix Q^ M^ in (19) only 
needs to be computed once. At each time step or iteration, only values of the nonlinear function 
need to be calculated without involving any numerical quadratures. The total complexity in each 
time step (or iteration) is reduced to 0{g{2rndof) + "^^^T^dof) flops. The efficiency can be further 
improved by combining with the DEIM, which will be discussed in Section 6. 

4 Error Estimates 

In this section, we analyze the numerical error of the POD-FEIC model (17). Since the obvious 
distinction between the POD-FEIC (17) and the POD-FEM (14) hes in the special spatial dis- 
cretization of the nonlinear term, we will consider the semidiscrete model (continuous in time) 
and only focus on errors caused by the spatial discretization and the POD truncation. To provide 
the analysis, we proceed in three steps: We begin by gathering a few preliminary results that will 
be used; We then prove an error estimate for the L^ projection of u in Lemma 4.6; Finally, we 
establish the approximation error of (17) in Theorem 4.1. 

For clarity of notation, in the sequel, we will denote by C, a generic constant that does not 
depend on the mesh size h and the number of POD basis functions r in the ROM. 

4.1 Step 1: Preliminary Results 

To derive the error estimation, we first make a few necessary hypotheses and present some prelim- 
inary results. For the solution u and the nonlinear function N{u) of (7), we assume: 

(HI) The solution u belongs to C\0,T; H'^+^{n)p\H^{n)), 

(H2) N{u) belongs to C(0,r; i/^+i(f^)), 

(H3) N{u) is locally Lipschitz, that is, let M = ||'^||l°°(^) + 1? there exists L = L(M) such that 

\N{u) -N{v)\ <L\u-v\, ioT sllu.v G {-M,M). (20) 

Based on the FE method theory, we have the following interpolation error [11]: 

Lemma 4.1 For < 7 < m + 1 and 1 < p < 00, if v e C{Q) f] Yl W^^^{K), there exists a 

KeTh 
constant C independent of h such that 

lu, _r^7,ii < r h^^^~^ WiiW ,1 oa^ 



The snapshots in our considerations are composed of the FE solutions ix^(x,t), which solves 
the following approximation problem: To find u^{x,t) G V^, such that. 



,v''^+a(u\v'') + (^N{u''),v'')=(f,v''), V^;'^Gy^ (22) 



~dt' 

with u^{x,0) — Uq{x) e V^, \/x £ fl. The FE approximation theory for semiUnear parabohc 
equations is well-developed (see, e.g., Chapter 14 in [68] and reference therein). One can easily 
modify the proof of Theorem 14.1 in [68] and obtain the following error estimate for the FE solution. 
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Lemma 4.2 Letu^ andu he the solutions of (22) and (7) under the assumptions of (HI) and (H3), 
respectively. With appropriately chosen u^ in the FE approximation, we have , with C — C{u^T), 



\u^{t) - u{t)\\ + h\\u^{t) - u{t)\\i < Ch 



m+l 



forte[0,T]. 



(23) 



To estimate the interpolation of nonlinear terms in the FEIC, we utilize an auxiliary 'euclidean' 
norm introduced in [69] on C{ft): 



h = 



ridof 
i=l 



Xi 



(24) 



For any x ^ ^hi since the space is of finite dimensions, we have the equivalence between ||x|| and 
llxlU on the reference element. With a straightforward homogeneity argument (or scaling argument 
[11]), we have the following lemma: 



Lemma 4.3 There exist two strictly positive constants c\ and C2 independent of h such that 



cih^WxWh < 11x11 <C2h2\\x\\h, 



(25) 



for allx ^ Vh. 



Recall that the FE solutions u^{x,t) are used as snapshots, 1-L ^ L? \s considered in the POD 
method and <^j{x) is the j-th POD basis. Besides the POD projection error in L^(0,r, L^(ri)) 
given in Proposition 2.1, 



r||A-,*)-E(«'(-,0,^i(-))^i(-)fdt = 5]A,- 



(26) 



j=i j>r 

we have the projection error in L'^{0,T,H^(ft)) norm as follows. 
Lemma 4.4 ([64]) The POD projection error in H^ norm satisfies 



j=i jyr 



(27) 



For the POD approximation, we have the following POD inverse estimate [43]: 

Lemma 4.5 Let Mr G W^^ with [Mr]jk = {(pj,(fk) be the POD mass matrix, Sr G W^^ with 
[Sr]jk — i^Ajk + (V(/^j, V(/?/e) bc the POD stiffness matrix, and \\ • ||2 denote the matrix spectral 
norm. Then, for all v ^Vr, the following estimates hold. 



IM 



r\\2 



iDr 



2 \\V\\l, 



\\v\\ < 

\\v\\i < ^J\\Sr\\2\\M;r^\\2\\v\\. 
Since we choose T-L = L^ in the POD method, both \\Mr\\2 and ||M~-^||2 are one. 



(28) 
(29) 
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4.2 Step 2: L^ Projection Error 

Next, we define the Lp' projection of u^ w^, from L^ to V^ such that 

u-w^,v^) =0, ^v^ eVj". 



-2 



(30) 



We have the following estimation of the L^ projection error. 

Lemma 4.6 The L? projection of u, w^, satisfies the following error estimations: 



f 

Jo 



f 


u-w'^ 


Jo 




v(i 


. - w'^) 



\ j>r J 



J ' 



(31) 
(32) 



j>r 



where C ^C{u,T). 
Proof 



u — w^ 



u — w^^u — w^ 



(30) 



yv^ evl". 



r ^ ' r 



It indicates, by Cauchy-Schwartz inequality, for all u^ G Vj^, 



u — w^ 



< u-v': . 



(33) 



(34) 



Decomposing u — v^ — u — u^ -^ {u^ — v^) and choosing v^ — V^u^ — ^ [u^^ (pj^ (pj in (34), by the 
triangular inequality and Proposition 2.1, we have 



u — w^ 



dt < C 



(26) 



< c 



(i: 



u — u 



u — u 



dt 



pT ^ 2 \ 



dt + E^i 



(35) 



Considering the FE approximation error (23), we have the bound for the L^ projection error in 

L2(0,r;L2(r2)) as (31). 



By the triangular inequality and adding and subtracting VrVp' — ^ {u^, ipj) ipj, we have 



f 

Jo 



V iu — w. 



dt < 



+ ||V(n'*-P^n'* 



V {Fru'' - w^ 



(27),(29) 

< c 



c I (||v(„-„'- 

dt + y^ ||(/^j 111 Aj + 115^112 / 11'^ -^ 
7^ -'o 



dt 






w^ 



h\\2 



..h\\2 
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where we use \\w^ — V^u^W < \\u — u^\\ in the last inequality. Considering the FE approximation 
error estimation (23), we have the bound for the L^ projection error in L'^{O^T; H^{Q)) as (32). 
This proved the lemma. 

4.3 Step 3: Main Results 

Finally, we discuss the main theoretical results of this paper, which is about the approximation 
property of the POD-FEIC (17). We first estimate the difference between the L^ projection of u, 
w^^ and the approximation solution u^ on a certain time interval in Lemma 4.7, which is bounded 
by the L? projection error. The conclusion is then extended to the whole time interval through a 
continuity argument in Lemma 4.8. By using the triangular inequality, we get the final estimation 
of approximation error u — u^ in Theorem 4.1. The crucial component in the analysis is the 
interpolation error of the nonlinear term [N(u) —X^N{u^)^v^^^ which can be decomposed in two 
different ways: 



and 



N{u) - I^N{u^),v^^ = (^N{u) - I^N{u),v^^ + (l^N{u) - I^N(u^),v^^ , (36) 

(^N{u) - I^N{u^),v^^ = (^N{u) - N{u^),v^^ + (^N{u^) - I^N{u^),v^^ . (37) 

The first approach has been used in [69], in which the first term on the RHS of (36) is bounded 
by the standard FE interpolation error under the smoothness assumption of N(u), and the second 
term can be estimated with the help of the local Lipschitz continuity assumption of N{u) and an 
auxiliary 'euclidean' norm. The second approach has been used in [19, 46], in which the first term 
on the RHS of (37) also appears in the standard FE discretization, thus can be treated as usual. 
The second term can be estimated by the FE interpolation error of N{u^), whose accuracy relies 
on the regularity of N{u^). It is determined by a hypothesis on u^ and a stronger smoothness 
assumption of N{u) than that required in the first approach. In this paper, we will follow the first 
approach (36). 

let V — v^ vn (7) and subtract it from (17), we have the error equation oi e — v!^ — u d^s follows. 

et,v^) ^a(e,v^) + (x^N{u^^) - N{u),v^) =0, Vi;^ G 1//^. (38) 



Let (j)^ — v!^ — w^ and r] — u — w^^ we have the decomposition of error, e — (fi — r]. Based on the 
definition of L^ projection (30), we have {rft^v^) — 0. Therefore, the error equation (38) can be 
rewritten as: 

((/.Jf,,, v^) + a (ckl v^) = a (r?, v^) + (n{u) - X^N{u^^),v';:) , Vt;,^ e V,\ (39) 

Lemma 4.7 Under assumptions (H1)-(H3), let u he the solution of (7) and u^ be the solution of 
(17) with the initial condition u^{-^0) = Yl^j =i{^o i ^r j)^r j ^^ (1^)^ respectively. Assume that, for 
some ti with < ti <T, we have 

||</.^(t)|Uoo(n) < ^, for anO<t< h. (40) 

Then it follows that 

mti)f + J^'\\^<t>Ht)fdt<C{u,T) Ih'-^ + \\Srhh'"'+' + J2\Ml^A ■ (41) 
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Proof First note that one can certainly find a ti to make (40) true by adjusting the mesh size 
h and the number of POD basis functions r. We choose h to be small enough and r to be large 
enough to ensure that ||t^ — tt;^||L°°(^) < I? then together with (40), we have both w^^ and u^ locate 
on the interval (— M, M), where M is defined in (H3). This allows us to take advantage of the local 
Lipschitz condition of N(u) on the interval t G [0,ti]. Next we start estimations. 
Let v^ — (f)^ in (39), we have 

(,/.^_„ </>:!) + a (,/>^ </>:?) = a (r7, </>:?) + [n{u) - X^N(u), <P^) 

+ {x^N{u) - I^N{w^), 0j!) + (x^N{w^) - l''N{4), </.J?) . (42) 

By the continuity (9) and coercivity (10) of a(-, •), and Cauchy-Schwartz inequality, we have 

^^ll^^f + /3||V</>:?f <«||Vr?||||V0,^|| + ||iV(n)-X'^iV(n)||||0:?|| 

+ \\l''N{u)-l''N{w'^)\\\\<p'^\\ + \\l''N{w^)-I^N{u'^)\\\\<p'^\\. (43) 

For the first term on the RHS of (43), by the Youngs' inequality, we have 

«l|Vr?|| ||V<^,^|| < ^llVryf + ^\\V<p'^f. (44) 

For the second term on the RHS of (43), Young's inequality yields 

||iV(„) - Z'.A.,„)|| ll^fll < "^""-f^W"' + ^. (45) 

Considering the FE interpolation error and the assumption (H2), we have 

\\Niu) -l''Niu)\\ < C/i-+i|iV(n)|H^+i(o). (46) 

For the third term on the RHS of (43), Young's inequality gives 

lii''iv(„) - i'.iv(.,;)ii ii^fii < "^"^'"'-f^'"-'"' + ME. (47) 

Employing the Lemma 4.3, the local Lipschitz condition, and the triangular inequality, we have 

\\l''N{u)-I^N{w^)\\ ? C2h-2\\l''Niu)-l''Niw'^)\\h ^'^='^^ C2h^\N{u) - N{w';^)\\h 

< Lc2q^||X'*n-<|| 

< Lc2c];^ (\\u-l'^u\\ + \\u-w'i^\\) . (48) 

By the regularity assumption (HI) of solution u and the FE approximability (Lemma 4.1), we have 

\\u-l''u\\ < C/i"*+VlH™+i(n)- (49) 

Thus, we have, for the third term on the RHS of (43), 

\\l''N{u)-l''N{w'^)\\ < C(/i"^+^ + ||r?||) . (50) 
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For the forth term on the RHS of (43), we use similar arguments to those for the third term and 
get 

(20) d.. u u.. (25) 



r ^r I 

l\\J.h\ 



< Lc2h2\\w]^ - u^Wh < Lc2C-y\\w^-u\ 



= Lc2qlCI|. (51) 

Substituting (44), (45), (46), (47), (50), and (51) in (43), we obtain 

^||0,^f + /3||V(/),^f < ^llVr/f + C(/i2-+2 + hf)+C*||0,^f, (52) 

where C* = Lc2C^^ + 1. By Gronwall's lemma, on the interval [0,ti], we have 

U^r{tl)\? + r/3||V0,^fd.<||0,^(O)f 

^^ati I ' (^^llVr/f + C/i2^+2 ^ c\\r]\A ds. (53) 



Considering ti < T and the choice of the initial condition which indicates ||0|?(O)|| = 0, the above 
inequality yields 



v'(ti)f + /'/3||V0,^fd. < e^*^^ (^||Vr;f + C/^2-+2 + q|^||2^d., 



(32) 






where C = C{u^T) independent of ti. 



Lemma 4.8 Suppose the order of finite elements m >1 for d — 1 and m > 2 for d>2, respec- 
tively. With the same conditions as those in Lemma 4-^, we have 

Proof Assume t| is the largest value that makes (40) true. If tl ^ T, it must be 

1 

2 
However, by the inverse inequality, 



HtmL^iCi)^^- (56) 



U'^{tl)h^^n)<Ch--2U'^{tm<Ch-2 U- + ,/fs;^/,-+i+ /J^||^^.||2aJ . (57) 

Then, for tti > 1 if rf = 1, and for tti > 2 if d = 2, 3, one can always find a h small enough and a r 
large enough such that ||0^(^i)||l°°(^) < ^- This contradicts with the assumption (56). Therefore, 
tl = T, that is, the conclusion (41) is true on the whole time interval [0,T]. 
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Theorem 4.1 Under assumptions (H1)-(H3), let u he the solution of (7) and u^ be the solution 
of (17) with the initial condition i/^(-,0) = X]^=i('^0 5 ^r jO^r ? ^^ (1^)^ respectively. The order of 
finite elements m > 1 for d — 1, m > 2 for d > 2. There exist positive numbers ho and ro such 
that, for h < Hq and r > r^, we have 



(58) 



\\4iT) - uiT)f + / ||V«^(t) - V«(t)f dt <cih^"' + WSrh /i''"+' + Yl ll<^ill? ^J ] 

where C = C{u,T) is a positive constant independent of h and Xj. 

Proof The conclusion follows the triangular inequality, Lemma 4.8, and Lemma 4.6. 

Remark 4.1 In Theorem ^--l, the estimation for the gradient of errors, ||Vix^ — Vi/||7^2(o,T;L2(^)) ^s 
optimal with respect to (w.r.t.) both the spatial discretization (^ 0{h^)) and the POD truncation 

(^ 0{J^-^^ ll^jlli '^j)y'' ^y Holder^s inequality, it is easy to show the error in L'^{O^T; H^(Q)) 

has the same optimality. Although not proven here, the error in L^(0,T; L^(ri)) is also optimal 

w.r.t. the spatial discretization (^ 0{h^^^)) and the POD truncation (^ (^{\/J2j>r^j))- ^^ 

shall verify it numerically in the next section. 

Remark 4.2 Theorem ^.i also holds on the FEIC approximation of the POD-G-ROM for the 
problem (8). The proof follows a similar argument. 

5 Numerical Tests 

The goal of this section is twofold: first, the proposed method will be validated by two problems 
appeared in interdisciplinary research. Accuracy and efficiency of the new approach are to be 
tested; second, the theoretical result in Section 4 will be numerically verified by another example. 
We mainly compare the results of the POD-FEIC (17) with those of the POD-FEM (14). The 
errors in two different norms are measured: 



£o{u,v) 



\ ^5^ll^(-'^^)-^(''^^)P' £i{u,v)^. ^ ^^ 



M 



u{',ti) -v{-,ti)\\l, 



where Sq is a discrete approximation of the error in L^(0,r; L^(ri)), while £i is a discrete approx- 
imation of the error in L^(0,r; i7^(ri)). For a fair comparison, all numerical tests reported in this 
paper are implemented on a PC with a 2.6GHz Intel Core 17 processor. As a criterion for efficiency, 
the CPU time of simulations is considered, which is the (on-line) time elapsed for the integra- 
tion only, excluding the (off-line) time for generating basis functions, precomputing matrices, and 
calculating errors, etc. 

5.1 Validation 

For the validation purpose, we consider two examples: a one-dimensional (ID) FitzHugh-Nagumo 
(F-N) system, which possesses a cubic polynomial nonlinearity; and a two-dimensional (2D) Buckley- 
Leverett equation (BLE), which has a non-polynomial nonlinearity. 
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FitzHugh-Nagumo System We first consider the simplified Hodgkin-Huxley model used in 
[18], which is a F-N system. The model is a nonlinear PDF system and describes the activation 
and deactivation dynamics of a spiking neuron. The system reads: 



dv 1 1 

— filAv v(v — O.ljf 1 — v) ^ — w - 

Ot II jl 

dw . 

— bv -^ ^w 

at 



c 
--C, 



v{x,0) = 0, w{x,0) = 0, 
v^{0,t) ^ -io{t), v^{L,t)^0, 



xe[0,L],t€[0,T], 

xe[0,L],te[0,T], 

xe[0,L], 
te[0,r], 



(59) 



where v{x^t) and w{x^t) are voltage and recovery voltage, respectively. Let u = [v^w]'^^ the weak 
form of the original system has the form of (7): 



^,v) +Mia(u,v) + M2 (u,v) + (N,v) = (f,v), 



(60) 



where a(u,v) = /i(Vu,Vv), N = iT;(v - 0.1)(1 - v),0 



M'" 



Ml 












and 



M, 



r n 


1 - 




M 


[ -b 


7 . 



We choose the same parameters as those utilized in [18], that is, L = 1, T = 8, 



-15t 



II = 0.015, h = 0.5, 7 = 2, c = 0.05 and the stimulus zo(^) = 50000^*^6" 

To generate snapshots, we use linear finite elements on a uniform mesh for spatial discretization 
with mesh size h — 1/512, and the Crank-Nicolson scheme for time integration with time step 
At = 1 X 10~^. Totally, 801 snapshots are collected and used to compute POD basis in L? space. 
Since the exact solution is unknown, we regard the FF solution as the benchmark. 

The errors and simulation time of the POD-FFM and the POD-FFIC are hsted in Table 1 
when the same number of POD basis functions, r, is used in both ROMs. It is seen that the FFIC 
discretization achieves same accuracy as that of the FF discretization, however, reduces the CPU 
time by a factor of 150. A comparison of the limit cycle projected onto the v — w plane among 
the FF solution, the POD-FFM (r = 5) solution, and the POD-FFIC (r = 5) solution is shown in 
Figure 1, which illustrates the correct limit cycle of the original system has been captured by the 
POD-FFIC when only 5 POD basis functions are used. 

Table 1: Frrors and the CPU time of the POD-FFM (14) and the POD-FFIC (17). Note that the 
POD-FFIC keeps the same accuracy as the POD-FFM, but saves CPU time by over 150 times. 





POD-FEM 


POD-FEIC 


1 


£oK,u'') 


£:i«,n'*) 


CPU time 


£oK,u^) 


5i«,n'*) 


CPU time 


3 


3.64e-03 


9.96e-02 


2.09e+02 


3.64e-03 


9.96e-02 


1.20 


5 


6.61e-04 


2.27e-02 


2.11e+02 


6.60e-04 


2.27e-02 


1.23 


7 


1.20e-04 


4.91e-03 


2.12e+02 


1.20e-04 


4.91e-03 


1.25 


9 


2.03e-05 


1.25e-03 


2.12e+02 


2.10e-05 


1.25e-03 


1.27 
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Figure 1: The limit cycle on the v — w plane: the POD approximation (r = 5) with either the 
FEM discretization (POD-FEM) or the FEIC discretization (POD-FEIC) coincides with that of 
the finite element solution (FEM). 

Buckley-Leverett Equation We then consider the 2D BLE, which is usually used to describe 
two phase flow in porous media with a gravitation pull in x-direction [50]. 




du ^ dfi (u) 

ot ox 

u{x, 0) = 



df2iu) 



dy 
u{x,t) = 0, 



0, xen,te[o,T], 



(61) 



where fi(u) 



u^+{i-uy^ 



and f2{u) = fi{u)[l - 5(1 



XG[0,L], 

X edn.te [0,r], 

. The weak form of the original system 



has the form of (8) with a{u^v) = jiiyu^Vv) and N{u) = [— /i, — /2]^- 

In this test, we choose the same parameters as used in [50]: fi = 0.1, T = 0.5, Vt = [—1.5, 1.5] x 
[—1.5, 1.5]. To generate snapshots, we use linear finite elements on a uniform triangular mesh for 
spatial discretization with mesh size h — 1/64, and the Crank-Nicolson scheme for time integration 
with time step At = 1 x 10~^. 51 snapshots are collected and used to compute POD basis in L? 
space. Due to the lack of exact solution, we also consider the FE solution to be the benchmark. 

The POD-FEM and the POD-FEIC approximation erros are hsted in Table 2 when r POD 
basis functions used in both ROMs. It is seen that the POD-FEIC obtains the same accuracy as 
that of the POD-FEM, but, decreases the CPU time by 40 times. A comparison among the FE 
solution, the POD-FEM, and the POD-FEIC result at t = 0.2 is shown in Figure 2. 

Table 2: Errors and the CPU time of the POD-FEM (14) and the POD-FEIC (17). Note that the 
POD-FEIC keeps the same accuracy as the POD-FEM, but saves CPU time by over 40 times. 





POD-FEM 


POD-FEIC 


/ 


eoiK^W") 


Si{u^,u'') 


CPU time 


£:o«,n'*) 


£:i«,n'^) 


CPU time 


5 


4.44e-03 


7.63e-02 


4.86e+03 


4.43e-03 


7.62e-02 


116 


10 


3.58e-04 


1.03e-02 


4.92e+03 


3.70e-04 


1.04e-02 


118 


15 


3.34e-05 


1.34e-03 


4.95e+03 


1.04e-04 


2.29e-03 


119 
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FEM sol. at t= 0.2 



POD-FEM error at t= 0.2 



POD-FEIC error at t= 0.2 






Figure 2: The simulations at t = 0.2: the FEM solution (left), the error of POD-FEM with r = 5 
(middle), and the error of POD-FEIC with r = 5 (right). 



5.2 Verification 

In this subsection, we will verify the theoretical results obtained in Section 4 through a problem 
with a non-polynomial nonlinearity. It is governed by the following equations: 



du 
'dt 



Au + sin(ix) = 


= /, 


xen,te [0,T], 


u = 


= 0, 


xedn,t€ [0, T] 


u{x,0) = 


= 9, 


X e^, 



(62) 



.,>.A,<lxlO-7. 



where / is determined by substituting a designated exact solution into the LHS of (62) and g 
is the exact solution at t = 0. In the test, we consider the problem in ID with exact solution 
ix = 0.5 sin(7rx) (10 tanh(x — t) + 1) on the domain Q = [0, 1] during the time interval t G [0, 1] . The 
exact solution is also our benchmark when calculating errors. The weak formulation is of the form 
(7) with a{u^v) — (Vu,Vv) and N{u) = sin(i/). We investigate the convergence properties of the 
POD-FEIC solution w.r.t. mesh size h and the number of POD basis r, respectively. 

To check the approximation order of the POD-FEIC solution w.r.t. h, we collect the finite 
element solution of the original system with linear elements (m = 1) and quadratic elements 
(m=2) respectively. Backward-Euler method is used for the time integration with a small time 
step At = 1 X 10~^. The number of POD basis functions are chosen such that V •. 
In this way, the spatial discretization error dominates the whole approximation property. 

The errors in both Sq and £i norms are shown in Table 3. Linear regressions indicate that, for 
linear elements, the order of convergence is 1.97 in £q norm and 0.98 in £i norm; for quadratic 
elements, the error convergence order is 2.94 in Sq norm and 1.95 in £i norm. The approximation 
orders are close to the optimal values in Theorem 4.1. Since the error analysis is asymptotic, as 
h decreases, the order approaches the optimal values (order m ^ 1 in. £q norm and order m in £i 
norm) . 

To check the approximation order of the POD-FEIC solution w.r.t. r, we collect the finite 
element solution of the original system when linear elements (p = 1) and quadratic elements 
(p=2) are utilized for spatial discretization, respectively, and backward-Euler method for the time 
integration. The mesh size h = 1/64 and the time step At = 1 x 10~^ are fixed. The number of 
POD basis functions are chosen so that \/j2j>r^J decays by a factor of 2. 

The errors in both £o and £i norms when linear elements are used are shown in Table 4 and 
those for quadratic elements are listed in Table 5. Linear regressions indicate that the convergence 
order of error in £i norm w.r.t. \Ylij>r ll^j ll/fi ^j ^^ ^-96 for linear elements and 1.00 for quadratic 
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Table 3: The approximation rate of the POD-FEIC w.r.t. h 



h 


m = 1 


m = 2 


fo«,«) £iK,u) 


fo«,M) fl«,M) 


1/8 
1/16 
1/32 
1/64 


2.09e-02 5.41e-01 
5.50e-03 2.81e-01 
1.39e-03 1.42e-01 
3.50e-04 7.12e-02 


2.41e-3 1.26e-l 
3.30e-4 3.42e-2 
4.22e-5 8.76e-3 
5.31e-6 2.20e-3 


order 


j^l.97 ^U.98 


j^2M ;^1.95 



elements. Wherein, the approximation order for quadratic elements is linear, which is optimal 
as indicated in Theorem 4.1. That of linear elements is close to 1 but slightly smaller than 1, 
which is because the discretization error tends to dominate the error as r incr eases. At the same 
time, note that the approximation order of error in £q norm w.r.t. A/^jy^^j is 1 for both linear 
and quadratic elements. Although not proven theoretically, the error in £q norm also converges 
optimally. 

Table 4: The approximation rate of the POD-FEIC w.r.t. r when linear finite elements are used. 









m = 


1 






y Ei>r ^i 


£o{u'},u) 


^Jj:j>rM\Hl^j 


£iiu';:,u) 


3 


4.13e-02 


4.60e-02 




6.26e-01 


6.21e-01 


4 


2.45e-02 


2.74e-02 




4.48e-01 


4.47e-01 


6 


9.07e-03 


1.02e-02 




2.23e-01 


2.31e-01 


7 


5.59e-03 


6.29e-03 




1.55e-01 


1.69e-01 


order 


- 


1.00 




- 


0.96 



Table 5: The approximation rate of the POD-FEIC w.r.t. r when quadratic finite elements are 
used. 









m = 


2 






V Ei>r- ^i 


£o{u'f,u) 


^/EjyrMl H^^j 


£'i(nj^,«) 


3 


4.15e-02 


4.60e-02 




6.28e-01 


6.19e-01 


4 


2.46e-02 


2.74e-02 




4.50e-01 


4.43e-01 


6 


9.18e-03 


1.02e-02 




2.25e-01 


2.21e-01 


7 


5.68e-03 


6.29e-03 




1.57e-01 


1.54e-01 


9 


2.20e-03 


2.41e-03 




7.44e-02 


7.36e-02 


order 


- 


1.00 




- 


1.00 



6 The Combination with the DEIM 

The discrete empirical interpolation method has been successfully applied in many nonlinear ROMs 
to reduce the computation complexity of the nonlinear terms [17, 18, 16, 34, 41, 67]. For a general 
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nonlinear function N{u{x^t)), the DEIM employs the ansatz 

r 

N{u) = J2^j{x)cj{t), (63) 

where tpj (x) is the j-th nonlinear POD basis obtained by applying the POD method on the nonlinear 
snapshots and r is the rank of the nonlinear POD basis. Based on the nonlinear POD basis vectors 
* = ['01, . . . , '0f], the DEIM optimally selects a set of interpolation points p := [pi, ... , pp] and 
approximates the nonlinear function by 

N{u) ^ *(PT*)-ipTN(i/), (64) 

where P is the matrix for selecting the corresponding p indices pi, . . . , pp. 

However, when the FE method is used for a spatial discretization, the nonlinear snapshot 
becomes {N{u{x,t)), h) as in the weak formulation (22). Therefore, generating nonlinear snapshots 
requires the inner product to be fulfilled, which costs lots of off-line time. In cases such as complex 
flows are studied, many interpolation points might be necessary to obtain a good approximation 
of the nonlinear terms. Since on-line simulations also need to evaluate the inner product over 
the elements sharing the selected DEIM points, the on-line time increases. These issues represent 
the main computational hurdles for applying DEIM in the FE setting, which, however, can be 
easily overcome by the POD-FEIC approach we proposed in Section 3. Indeed, in the FEIC, the 
nonlinear snapshot can be chosen to be the value vector of the nonlinear function, A^(Qa(t)), 
thus, no any evaluation of inner product needed. Replacing the nonlinear function with the DEIM 
approximation (64) in the POD-FEIC (17), we get the POD-FEIC-DEIM model, which has a more 
efficient approximation of the nonlinear term 

AfpEic-DEiM = QTM'**(PT*)-ipTN(Qa(t)). (65) 

It is seen that, in on-line simulations, one only need to calculate the nonlinear functions at 
p selected DEIM points, which doesn't involve any numerical quadratures. Therefore, the POD- 
FEIC-DEIM improves the computational efficiency over the POD-FEIC, which also outperforms 
the POD-FEM. We will demonstrate the effectiveness of the new approach by considering the first 
two examples in Section 5 again. 

FitzHugh-Nagumo System We revisit the ID F-N model used in Section 5.1. Totally, 801 
nonlinear snapshots of N(v) = -v{v — 0.1) (1 — v) are generated and used in the DEIM. The 
nonlinear function is then approximated by the DEIM basing on p selected interpolation points. In 
this test, we consider the POD-FEIC-DEIM generated by the first r = 5 POD basis functions and 
investigate the numerical performance of the new model by varying the number of interpolation 
points. 

Errors in Sq norm and the CPU time for simulations are listed in Table 6. As p increases, the 
POD-FEIC-DEIM result approaches to that of the POD-FEIC {£o{u^,u^) = 6.60 x 10"^ when 
r = 5). It is also seen from Table 6 that when j> = 3, the POD-FEIC-DEIM error is 9 times larger 
than that of the POD-FEIC. However, as the limit cycle on the v — w plane shown in Figure 3, 
the difference mainly occurs at the beginning of the simulation, from t = to t = 1. After the 
transient interval, the limit cycles of the POD-FEIC and the POD-FEIC-DEIM coincide with each 
other. What's more, the CPU time of the POD-FEIC-DEIM is lower than that of the POD-FEIC, 
although not significantly due to the small size of the tested problem. 

19 



Nonlinear Model Reduction for Semilinear Parabolic Equations 



Z. Wang 



Table 6: Errors and the CPU time of the POD-FEIC-DEIM model with r 
interpolation points for nonlinear function approximation. 



5 POD basis and p 



P £o{u 



J~^) 



CPU time 



3 5.94e-03 

5 3.00e-03 

7 9.66e-04 

9 6.63e-04 



1.02 
1.05 
1.09 
1.09 



0.2^ 
0.18 
0.16 
0.14 
0.12 
^ 0.1 
0.08 
0.06 
0.04 
0.02 

-8.5 



- POD-FEIC 
- - POD-FEIC-DEIM p=3 




0.5 

V 



1.5 



Figure 3: The hmit cycle on the v - w plane: the POD-FEIC results (r = 5) and the POD-FEIC- 
DEIM approximation (r = 5 and p — ?>). 

Buckley-Leverett Equation We also consider the 2D BLE problem utilized in Section 5.1. 
Totally, 51 nonlinear snapshots of fi{u) and f2{u) are generated, and then used in the DEIM 
for selecting p interpolation points, respectively. We also consider the POD-FEIC-DEIM model 
generated by the first r = 5 POD basis functions and investigate its numerical behavior by varying 
p. 

Errors in £q norm and the CPU time for simulations are listed in Table 7. Note that for the 
POD-FEIC when r = 5, the error £^{u^,u^) = 4.43 x 10"^ and CPU time equals 116 5, the POD- 
FEIC-DEIM model achieves a close accuracy by only using p — ^ interpolation points, which also 
improves the computational efficiency of the POD-FEIC by 4 times. Figure 4 shows the distribution 
of 20 firstly selected DEIM interpolation points (left), and the difference between the POD-FEIC 
and the POD-FEIC-DEIM when p = 5 (middle) and j9 = 20 (right). 



Table 7: Errors and the CPU time of the POD-FEIC-DEIM model with r 
interpolation points. 



5 POD basis and p 



p 


^o«,n'^) 


CPU time 


5 


4.81e-03 


27.79 


10 


4.64e-03 


27.89 


15 


4.59e-03 


28.00 


20 


4.54e-03 


28.57 
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DEIM interpolation points 



POD-FEIC-DEIM (p=5) error at t= 0.2 



is -1 -0.! 



5 0.5 1 1.5 




POD-FEIC-DEIM (p=20) error at t= 0.2 




Figure 4: Distribution of DEIM interpolation points for fi{u) and f2{u) (left), the difference 
between POD-FEIC and POD-FEIC-DEIM at t = 0.2 when p = 5 (middle) and p = 20 (right). 

7 Conclusions 

As a first step of our investigations on efficient finite element discretization algorithms for nonlinear 
model reduction techniques, we develop the finite element method with interpolated coefficients for 
nonlinear POD-ROMs. Comparing with the standard finite element discretization, the proposed 
approach is computationally more efficient because the coefficient matrices in the nonlinear terms 
can be precomputed and there is no any numerical quadratures needed during the simulation. 
The proposed method also achieves the same accuracy as that of the standard finite element 
discretization when nonlinear functions satisfy certain smoothness assumptions. The proposed 
approach is a more suitable base for the discrete empirical interpolation method. Combining the 
FEIC with the DEIM will further reduce the computational complexity for evaluating nonlinear 
terms. 

We plan to continue investigating several research avenues: we will first extend the proposed 
approach to nonlinear closure models we developed for complex flows in [71, 72]. Some realistic 
engineering application problems will be tested. Second, we will develop a rigorous error analysis of 
the model reduction approach that combining the FEIC and the DEIM. Finally, we plan to design 
an adaptive algorithm for the proposed method, which will be used to supervise the approximation 
error control. 
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